This document will explore the distribution of MLGs across years as well as filtering strategies.
library(PramCurry)
## Loading required package: poppr
## Loading required package: adegenet
## Loading required package: ade4
## ==========================
## adegenet 1.4-2 is loaded
## ==========================
##
## - to start, type '?adegenet'
## - to browse adegenet website, type 'adegenetWeb()'
## - to post questions/comments: adegenet-forum@lists.r-forge.r-project.org
##
##
## This is poppr version 1.1.2.99.70. To get started, type package?poppr
library(reshape2)
library(ggplot2)
library(ape)
library(poppr)
library(adegenet)
library(igraph)
##
## Attaching package: 'igraph'
##
## The following object is masked from 'package:ape':
##
## edges
library(animation)
options(stringsAsFactors = FALSE)
data(ramdat)
data(pop_data)
data(myPal)
sessionInfo()
## R version 3.1.2 (2014-10-31)
## Platform: x86_64-apple-darwin13.4.0 (64-bit)
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] animation_2.3 igraph_0.7.1 ape_3.1-4.5 ggplot2_1.0.0
## [5] reshape2_1.4 PramCurry_1.0 poppr_1.1.2.99-70 adegenet_1.4-2
## [9] ade4_1.6-2
##
## loaded via a namespace (and not attached):
## [1] assertthat_0.1 bitops_1.0-6 boot_1.3-11
## [4] caTools_1.17.1 colorspace_1.2-4 digest_0.6.4
## [7] dplyr_0.2 evaluate_0.5.5 fastmatch_1.0-4
## [10] formatR_1.0 gdata_2.13.3 ggmap_2.3
## [13] grid_3.1.2 gtable_0.1.2 gtools_3.4.1
## [16] htmltools_0.2.6 httpuv_1.3.0 knitr_1.6
## [19] labtools_0.4.1 lattice_0.20-29 lubridate_1.3.3
## [22] mapproj_1.2-2 maps_2.3-9 MASS_7.3-34
## [25] Matrix_1.1-4 memoise_0.2.1 munsell_0.4.2
## [28] nlme_3.1-117 parallel_3.1.2 pegas_0.6
## [31] permute_0.8-3 phangorn_1.99-7 plotrix_3.5-7
## [34] plyr_1.8.1 png_0.1-7 proto_0.3-10
## [37] Rcpp_0.11.2 RgoogleMaps_1.2.0.6 rjson_0.2.14
## [40] RJSONIO_1.3-0 rmarkdown_0.3.10 scales_0.2.4
## [43] shiny_0.10.1 stringr_0.6.2 tools_3.1.2
## [46] vegan_2.0-10 xtable_1.7-4 yaml_2.1.13
First thing to do is to create the full MSN across years. Note that the get_layout function will preserve the layout (“MASTER”) so that it can be subsetted.
setpop(ramdat) <- ~ZONE2
options(digits = 10)
(newReps <- other(ramdat)$REPLEN)
## PrMS6A1 Pr9C3A1 PrMS39A1 PrMS45A1 PrMS43A1
## 3 2 2 4 4
newReps[3] <- 4
(newReps <- fix_replen(ramdat, newReps))
## PrMS6A1 Pr9C3A1 PrMS39A1 PrMS45A1 PrMS43A1
## 3.00000 2.00000 4.00001 4.00000 4.00000
b.msn.full <- bruvo.msn(ramdat, replen = newReps, showplot = FALSE,
include.ties = TRUE)
goodSeed <- 19
goodSeed2 <- 17
goodSeedAlt <- 14
# for (i in goodSeed:100){
set.seed(goodSeed)
# set.seed(i)
MASTER <- get_layout(b.msn.full$graph)
plot_poppr_msn(ramdat, b.msn.full, gad = 10, palette = funky, mlg = TRUE,
layfun = MASTER, nodebase = 1.75, vertex.label.font = 2,
quantiles = FALSE,
vertex.label.color = "firebrick", inds = "none")
# prompt <- paste("save", i, "as seed?")
# accept <- readline(prompt)
# if (substr(accept, 1, 1) == "y"){
# theSeed <- i
# break
# }
# }
Now we will see how the MSNs are plotted per year. For this, we will do two things:
saveGIF({
plot_poppr_msn(ramdat, b.msn.full, gad = 10, palette = funky, mlg = TRUE,
layfun = MASTER, nodebase = 1.75, vertex.label.font = 2,
quantiles = FALSE,
vertex.label.color = "firebrick")
for (i in levels(pop(ramdat))){
population_subgraph(ramdat, b.msn.full,
quantiles = FALSE,
gad = 10,
palette = funky,
mlg = TRUE,
layfun = MASTER,
nodebase = 1.75,
vertex.label.font = 2,
cutoff = 0.148,
nodelab = 1000,
vertex.label.color = "firebrick",
sublist = i)
}
}, movie.name = "by_year.gif", interval = 1, ani.height = 1000, ani.width = 1000)
## Executing:
## 'convert' -loop 0 -delay 100 Rplot1.png Rplot2.png Rplot3.png
## Rplot4.png Rplot5.png Rplot6.png Rplot7.png Rplot8.png
## 'by_year.gif'
## Output at: by_year.gif
## [1] FALSE
# pdf("by_year.pdf", width = 10, height = 10)
plot_poppr_msn(ramdat, b.msn.full, gad = 10, palette = funky, mlg = TRUE,
layfun = MASTER, nodebase = 1.75, vertex.label.font = 2,
quantiles = FALSE,
vertex.label.color = "firebrick")
for (i in levels(pop(ramdat))){
population_subgraph(ramdat, b.msn.full,
quantiles = FALSE,
gad = 10,
palette = funky,
mlg = TRUE,
layfun = MASTER,
nodebase = 1.75,
vertex.label.font = 2,
cutoff = 0.148,
nodelab = 1000,
vertex.label.color = "firebrick",
sublist = i)
}
# dev.off()
# Just to see once more with only allowing 1 mutational step for Bruvo's distance
plot_poppr_msn(ramdat, b.msn.full,
quantiles = FALSE,
gad = 10,
palette = funky,
mlg = TRUE,
layfun = MASTER,
nodebase = 1.75,
vertex.label.font = 2,
cutoff = 0.051,
vertex.label.color = "firebrick")
setpop(ramdat) <- ~ZONE1
# zonePal <- get_year_pal(ramdat)
b.msn.zone1 <- bruvo.msn(ramdat, replen = newReps, showplot = FALSE,
include.ties = TRUE)
plot_poppr_msn(ramdat,
b.msn.zone1,
gadj = 10,
palette = funky,
mlg = TRUE,
layfun = MASTER,
nodebase = 1.75,
vertex.label.font = 2,
quantiles = FALSE,
vertex.label.color = "firebrick")
for (i in levels(pop(ramdat))){
population_subgraph(ramdat,
b.msn.zone1,
quantiles = FALSE,
gad = 10,
palette = funky,
mlg = TRUE,
layfun = MASTER,
nodebase = 1.75,
vertex.label.font = 2,
cutoff = 0.148,
vertex.label.color = "firebrick",
sublist = i,
nodelab = 1000,
inds = "none")
}
setpop(ramdat) <- ~ZONE2
# zonePal <- get_year_pal(ramdat)
b.msn.zone <- bruvo.msn(ramdat, replen = newReps, showplot = FALSE,
include.ties = TRUE)
zdegs <- igraph::degree(b.msn.zone$graph)
names(zdegs) <- igraph::V(b.msn.zone$graph)$label
sort(zdegs)
## MLG.26 MLG.19 MLG.65 MLG.64 MLG.55 MLG.58 MLG.43 MLG.7 MLG.40 MLG.11
## 1 1 1 1 1 1 1 1 1 1
## MLG.42 MLG.38 MLG.9 MLG.62 MLG.10 MLG.37 MLG.31 MLG.61 MLG.1 MLG.63
## 1 1 1 1 1 1 1 1 1 1
## MLG.20 MLG.69 MLG.67 MLG.70 MLG.32 MLG.25 MLG.5 MLG.44 MLG.36 MLG.35
## 2 2 2 2 2 2 2 2 2 2
## MLG.13 MLG.2 MLG.45 MLG.8 MLG.27 MLG.16 MLG.41 MLG.3 MLG.51 MLG.66
## 2 2 2 2 3 3 3 3 3 3
## MLG.33 MLG.59 MLG.6 MLG.57 MLG.56 MLG.50 MLG.34 MLG.39 MLG.4 MLG.24
## 3 3 3 3 3 3 3 3 4 4
## MLG.28 MLG.68 MLG.60 MLG.52 MLG.46 MLG.54 MLG.15 MLG.47 MLG.21 MLG.29
## 4 4 4 4 4 4 4 4 5 5
## MLG.14 MLG.30 MLG.18 MLG.12 MLG.49 MLG.23 MLG.53 MLG.48 MLG.22 MLG.17
## 5 5 5 5 5 6 6 6 7 8
edges_to_remove <- E(b.msn.zone$graph)[E(b.msn.zone$graph)$weight > 0.05]
clusts <- clusters(delete.edges(b.msn.zone$graph, edges_to_remove))
(names(clusts$membership) <- V(b.msn.zone$graph)$label)
## [1] "MLG.22" "MLG.27" "MLG.26" "MLG.4" "MLG.24" "MLG.21" "MLG.16"
## [8] "MLG.20" "MLG.19" "MLG.41" "MLG.17" "MLG.28" "MLG.29" "MLG.23"
## [15] "MLG.68" "MLG.69" "MLG.3" "MLG.51" "MLG.14" "MLG.65" "MLG.64"
## [22] "MLG.67" "MLG.70" "MLG.60" "MLG.52" "MLG.66" "MLG.30" "MLG.46"
## [29] "MLG.33" "MLG.32" "MLG.59" "MLG.55" "MLG.58" "MLG.53" "MLG.18"
## [36] "MLG.54" "MLG.43" "MLG.12" "MLG.48" "MLG.7" "MLG.25" "MLG.5"
## [43] "MLG.40" "MLG.11" "MLG.6" "MLG.44" "MLG.42" "MLG.38" "MLG.9"
## [50] "MLG.62" "MLG.36" "MLG.15" "MLG.10" "MLG.37" "MLG.35" "MLG.31"
## [57] "MLG.13" "MLG.57" "MLG.61" "MLG.47" "MLG.1" "MLG.49" "MLG.56"
## [64] "MLG.2" "MLG.50" "MLG.63" "MLG.34" "MLG.45" "MLG.39" "MLG.8"
clustPal <- clusts$csize
theClusts <- clustPal > 3
clustOrder <- order(clustPal, decreasing = TRUE)
clustPal[clustOrder][theClusts[clustOrder]] <- RColorBrewer::brewer.pal(sum(theClusts), "Set1")
clustPal[clustOrder][!theClusts[clustOrder]] <- gray.colors(sum(!theClusts))
nodeList <- lapply(1:length(clustPal), function(x) which(clusts$membership == x))
# saveGIF({
plot_poppr_msn(ramdat, b.msn.zone, gad = 10, palette = funky, mlg = TRUE,
layfun = MASTER, nodebase = 1.75, vertex.label.font = 2,
quantiles = FALSE,
vertex.label.color = "firebrick",
mark.groups = nodeList[theClusts],
mark.border = clustPal[theClusts],
mark.col = transp(clustPal[theClusts], 0.05),
mark.expand = 2,
mark.shape = 0)
for (i in levels(pop(ramdat))){
population_subgraph(ramdat,
b.msn.zone,
gadj = 10,
palette = funky,
mlg = TRUE,
layfun = MASTER,
nodebase = 1.75,
vertex.label.font = 2,
quantiles = FALSE,
vertex.label.color = "firebrick",
sublist = i,
inds = "none",
nodelab = 1000)
}
# }, movie.name = "by_zone2.gif", interval = 1, ani.height = 1000, ani.width = 1000)
mainMLGs <- order(table(ramdat@mlg), decreasing = TRUE)[1:10]
pdf("by_zone2.pdf", width = 352/2.54/10, height = 352/2.54/10, pointsize = 40)
plot_poppr_msn(ramdat, b.msn.zone, gad = 10, palette = funky, mlg = TRUE,
layfun = MASTER, nodebase = 1.75, vertex.label.font = 2,
quantiles = FALSE,
# inds = mainMLGs, nodelab = 1000,
vertex.label.color = "firebrick",
mark.groups = nodeList[theClusts],
mark.border = clustPal[theClusts],
mark.col = transp(clustPal[theClusts], 0.05),
mark.expand = 2,
mark.shape = 0)
for (i in levels(pop(ramdat))){
population_subgraph(ramdat,
b.msn.zone,
gadj = 10,
palette = funky,
mlg = TRUE,
layfun = MASTER,
nodebase = 1.75,
vertex.label.font = 2,
quantiles = FALSE,
vertex.label.color = "firebrick",
sublist = i,
inds = "none",
nodelab = 1000)
}
dev.off()
## pdf
## 2
zone2pal <- char2pal(setpop(ramdat, ~ZONE2)@pop.names, funky)
setpop(ramdat) <- ~ZONE2/Pop
zone2yearpal <- zone2pal[sub("_.+?$", "", ramdat@pop.names)]
names(zone2yearpal) <- ramdat@pop.names
# zonePal <- get_year_pal(ramdat)
b.msn.zone.year <- bruvo.msn(ramdat, replen = newReps, showplot = FALSE,
include.ties = TRUE)
# saveGIF({
# pdf("by_zone2year.pdf", width = 10, height = 10)
plot_poppr_msn(ramdat, b.msn.zone.year, gad = 10, palette = funky, mlg = TRUE,
layfun = MASTER, nodebase = 1.75, vertex.label.font = 2,
quantiles = FALSE,
vertex.label.color = "firebrick")
for (i in sort(levels(pop(ramdat)))){
population_subgraph(ramdat,
b.msn.zone.year,
gadj = 10,
mlg = TRUE,
layfun = MASTER,
nodebase = 1.75,
vertex.label.font = 2,
quantiles = FALSE,
palette = funky,
vertex.label.color = "firebrick",
sublist = i,
inds = "none",
nodelab = 1000,
main = paste("\n\n", i))
}
# dev.off()
# }, movie.name = "by_zone2year.gif",
# interval = 1, ani.height = 1000, ani.width = 1000)